function [R rr rm rs rrsig] = re (y, x, xm)

   ## usage:  [R rr rm rs rrsig] = re (y, x, xm)
   ## 
   ## reduction of error for y simulating x (relative to xm)
   
   ## Author:  GB
   ## Description:  reduction of error

   global alpha
   if isempty(alpha), alpha = 0.05 ; endif

   pkg load nan

   [nx, mx] = size(x) ;
   [ny, my] = size(y) ;

   if (nx != ny)
      error("x and y row mismatch") ;
   endif 
   if (mx == 1)
      x = repmat(x, 1, my) ;
   endif 

   I = !isnan(d = y-x) ; n = sum(I) ; d(!I) = 0 ;
   if !any(I), R = rr = rs = rm = p = ci1 = ci2 = nan ; return ; endif

   ve = sumsq(d) ./ n ;

   if (nargin > 2)
      xx = x - xm ;
   else
      xx = center(x) ;
   endif
   I = !isnan(xx) ;
   vx = sumsq(xx) ./ n ;

   wrn = warning ;
   warning("off", "Octave:divide-by-zero") ;

   R = 100 * (1 - ve./vx) ;

   if (nargout == 1)
      warning(wrn) ;
      return ;
   else
      # xm = mean(x) ; ym = mean(y) ;
      # x = center(x) ; y = center(y) ;
      rs = std(y)/std(x) ;
      rm = -100 * ((mean(y)-mean(x))/std(x))^2 ;
      # rr = 100 * (2*corrcoef(y,x) - rs) * rs ;
      rr = 100 * corrcoef(y, x) ;
      d = n - 2 ; tc = tinv(1 - alpha, d) ; rrsig = tc ./ sqrt(d + tc.^2) ;
      if (abs(imag(rr)) < eps), rr = real(rr) ; endif 
   endif 
   warning(wrn) ;
   return ;

   if (ve <= vx)
      R = 100 * (1 - ve./vx) ;
   else
      R = 100 * (vx./ve - 1) ;
   endif 

endfunction
